In this appendix, we review some familiar topics in numerical linear algebra that are used in this text and introduce some topics that are not common knowledge. The material here is not necessary to understand the book (assuming the reader will accept certain algorithms, specifically intersection multiplicity, without knowing the details). It is assumed the reader of this appendix will have some prior knowledge of linear algebra and matrices. The expert in numerical linear algebra via Mathematica can jump to A.3. The recommended source for further reading is the wonderful linear algebra text by Gilbert Strang (1976).

As an aside, in my younger days, I was a very abstract, pure-mathematics-only person. My interest in numerical mathematics originated in teaching linear algebra from Strang’s textbook.

A.1 Basic Real Linear Algebra

Like classical algebraic geometry, classical linear algebra lives in the complex domain. A main tool is the eigenvalue-eigenvector decomposition, which requires complex numbers—even for simple real matrices. However, modern real numerical linear algebra is based on the singular value decomposition, which gives real singular values even for complex matrices. First, we review Mathematica’s linear algebra implementation.

A real numerical matrix is a m×n array, which is a list of m lists of n machine numbers.

A=RandomReal[{-9,9},{3,4}]

A//MatrixForm

Columns and rows can be extracted. To do so, I prefer Part or Extract rather than Take.

c2= Extract[A,{All,2}]

A[[All,2]]

Take[A,All,{2}]//MatrixForm

Extract[A,2]//MatrixForm

r2=A[[2]]

Extract[A,2]

While the proper notation for a column is a list of one-element lists (as above), for multiplication purposes, the Wolfram Language will interpret a list of n numbers as either a row or column, depending on the context.

c2.A

Next, we review Strang’s four fundamental subspaces (1976) associated with a real m×n matrix A:

The “row space” is the the subspace of (ℝ^n) generated by the rows of A (that is, the set of all linear combinations of rows of A).

The “column space” is the subspace of (ℝ^m) generated by all columns of A.

The “null space” is the subspace of (ℝ^n) generated by all vectors v with A.v=0.

The “left null space” is the subspace of (ℝ^m) generated by all vectors w with w.A=0.

Two real vectors of the same length (v and w) are called “orthogonal” if v.w=0. I emphasize that this definition only works in the real case, not in the complex case. The orthogonal complement of a subspace V of (ℝ^n) is the subspace of (ℝ^n) consisting of all vectors w, which satisfy v.w=0 for all v in V. The orthogonal complement of V is denoted (V^⊥). Note ((((V^⊥)))^⊥)=V.

Strang’s fundamental theorem of real linear algebra is, “The orthogonal complement of the row space is the null space and conversely. The orthogonal complement of the column space is the left null space and conversely” (1976).

Again, this only works for real linear algebra, not in the complex setting (similarly for the next definition).

A square real matrix A is orthogonal if all rows (and equivalently all columns) are orthonormal—that is, v.w=0 when v and w are different rows (columns) and v.v=1 for all rows (columns). Alternatively, A is orthogonal when the inverse matrix to A is just the transposition of A, (A^(-1))=(A^T).

R= RotationMatrix[RandomReal[{0,2Pi}],RandomReal[{-1,1},3]]

Transpose[R]

Inverse[R]

A.2 The Singular Value Decomposition

The key idea in real numerical linear algebra is the singular value decomposition (SVD), which gives real results for real matrices (as opposed to the eigenvalue-eigenvector decomposition, which often gives complex results even for real matrices).

Given a real m×n matrix A, there exists an orthogonal m×m matrix U, a diagonal m×n matrix S with decreasing non-negative entries and an orthogonal n×n matrix U with A=U.S.(V^⊤).

A=RandomReal[{-4,4},{3,4}];
A//MatrixForm

{U,S,V}=SingularValueDecomposition[A];
U //MatrixForm

S //MatrixForm

V //MatrixForm

Chop[U.Transpose[U]]

Chop[V.Transpose[V]]

U.S.Transpose[V]//MatrixForm

The entries of S are called the singular values. One can do this with exact matrices, and if A is square and invertible, all singular values will be nonzero. If A is not invertible, some of the singular values will be zero. In fact, for an exact matrix, the number of nonzero singular values is the dimension of the row space (which is equal to the dimension of the column space). This number is called the “rank”. Further, for any m×n exact matrix A, the dimension of the null space is n-r and the dimension of the left null space is m-r.

For numerical matrices, there may be difficulties in determining the dimensions of the subspaces. The SVD allows us to resolve this problem. We pick a tolerance ϵ>0 (which is generally quite small but larger than the machine precision) and declare that any singular value s<ϵ is equal to zero. The number of singular values greater than ϵ is the rank r. The null space is then generated by the last n-r columns of V while the left null space is generated by the last m-r columns of U. The row space is generated by the first r columns of V while the column space is generated by the first r columns of U.

As an example, we set ϵ=(10^(-10)), which is actually a larger-than-normal tolerance.

B={{-1.2125910411564504`,-2.843903596715471`,-2.618517514835883`,-2.7219449070534747`},{-1.001767330326654`,-1.4814251186573628`,-3.266788014281124`,3.384583512352174`},{0.8773941003711174`,2.2586088618735145`,1.6393429153186834`,3.272955439991057`}};
B//MatrixForm

{U,S,V}=SingularValueDecomposition[B];
U//MatrixForm

S //MatrixForm

V//MatrixForm

S[[3,3]]=0;
S //MatrixForm

U.S.Transpose[V] //MatrixForm

The display version of this matrix (with the small singular value set to 0) does not vary from the original B, but the rank of this matrix is 2. The null space of this matrix is two dimensional and generated by the last two columns of V, as is the left null space, which is generated by the last column of U. The column space is generated by the first two columns of U, and the row space is generated by the first two rows of (V^T).

NS=Take[V,All,-2]; 
NS //MatrixForm
err1=B.NS

u0=Extract[U,{All,3}]

err2=u0.B

CS=Take[U,All,2];
CS//MatrixForm

Take[Transpose[V],2]//MatrixForm

One standard application is the pseudoinverse. Given A=U.S.(V^⊤), the pseudoinverse is obtained by killing diagonal entries in S less than the tolerance, transposing and replacing the remaining nonzero values by their reciprocals (call this ST). Then the pseudoinverse is P=V.ST.(U^⊤).

P=PseudoInverse[B,Tolerance->10^-10];
P//MatrixForm

B1=B.P;
B1//MatrixForm

B1.B1

B1.CS//MatrixForm

This is used in least squares computations. For example, the Gauss–Newton method refers to a method for finding least squares solutions of over-determined systems.

A special case with three equations in two variables is the following modification of the two-variable version of the Gauss–Newton method, replacing the matrix Inverse with the PseudoInverse.

gaussNewton3[f_,g_,h_,p0_,x_,y_]:=Module[{Q,Q0,p1,p2,J0,n},
n=30;
p1=p0;
Catch[Do[Q={f,g,h}/.Thread[{x,y}->p1];
J0=D[{f,g,h},{{x,y}}]/.Thread[{x,y}->p1];
p2=p1-PseudoInverse[J0].Q;
If[Norm[p1-p2]<1.*^-7,Throw[p2]];
p1=p2,{n}]];
p2]

Suppose we want to find a singular point of g=0 using the classical definition that a singular point is where g, D[g,x] and D[g,y] all vanish, starting from an approximation taken from a plot.

g=-0.22495692374518167`-1.1013674560719149` x+8.256356096694816` (x^2)-6.3764001275126425` (x^3)+1.3683563251360158` y-4.696428009684947` x y+7.522388070354437` (x^2) y-0.0035046560139825322` (y^2)-3.591849285015912` x (y^2)+0.7352364129554713` (y^3);
p={1.58,3.58};
{g,D[g,x],D[g,y]}/.Thread[{x,y}->p]

q=gaussNewton3[g,D[g,x],D[g,y],p,x,y]
{g,D[g,x],D[g,y]}/.Thread[{x,y}->q]

This is close enough to meeting the classical definition of a singular point! Much of the time, we will only be interested in the numerical rank of a matrix.

SingularValueList[B]

The Wolfram Language uses a default tolerance of (Subscript[s, 1])*100*(10^(-p)) , where p=$MachinePrecision and (Subscript[s, 1]) is the largest singular value and smaller singular values are listed as 0. This can be reset by the user, but if you are finding the rank and it is possible that a matrix is the zero matrix, you will get a list of zeros.

M2={{1,0,0,0},{1,1,0,0},{2,3,0,1.*^-13}};
Select[SingularValueList[M2],#>1.*^-12&]

Programming note: If M is an exact matrix, SingularValueList and SingularValueDecomposition will attempt an exact calculation of these quantities. For a large matrix (such as the matrices in the next two sections), this can take a large amount of time and not give useful results. Thus it is probably a good idea to use these in the form SingularValueList[N[M]] or SingularValueDecomposition[N[M]] unless you are sure you have a numerical matrix, that is, at least one entry contains a decimal point.

A.3 Sylvester Matrices

The material here works for polynomials of many variables. We will restrict it to the two-variable case to keep it simple.

This returns a list of rules (as it did in Chapter 9). However, perhaps a more convenient method would be to use Association.

For example, to find the total degree of a bivariate polynomial, we use:

tDeg[f_,x_,y_]:=Max[Table[Total[p],{p,Keys[Association[CoefficientRules[f,{x,y}]]]}]];

pExps[d_]:=Module[{pex},
pex[0]={{0,0}};
pex[n_]:= pex[n]=Join[pex[n-1],Table[{n-j,j},{j,0,n}]];
pex[d]]

pExps[3]

Given a bivariate polynomial f(x,y), we construct an association of coefficients by the following procedure. Since we often need the degree, we can save time by returning the degree as well as the association.

fAssoc[f_,x_,y_]:=Module[{A,B,d,FA},
A=Association[CoefficientRules[f,{x,y}]];
d=Max[Table[Total[p],{p,Keys[A]}]];
FA=Association[Table[If[MissingQ[A[p]],p->0,p->A[p]],{p,pExps[d]}]];
{FA,d}];

shiftFA[FA_,q_,n_]:=Module[{S,K,d},
K=Keys[FA];
S=Association[Table[p+q->FA[p],{p,K}]];
Association[Table[If[MissingQ[S[p]],p->0,p->S[p]],{p,pExps[n]}]]];

sylD1[f_,n_,x_,y_]:=Module[{FA,d},
{FA,d}=fAssoc[f,x,y];
If[d>n,Print["Must have n≥deg(f)"];Abort[]];Table[Values[shiftFA[FA,q,n]],{q,pExps[n-d]}]]

Finally, the general Sylvester matrix for two polynomials is given by the following (note that it would be easy to modify this procedure to handle an arbitrary number of polynomials, but that won’t be necessary for us):

sylvesterMatrix[f_,g_,n_,x_,y_]:=Join[sylD1[f,n,x,y],sylD1[g,n,x,y]];

f=x-y+x^4;
g=x-y+y^4;
S=sylvesterMatrix[f,g,5,x,y];
S//MatrixForm

One application of sylD1 is polynomial multiplication. Note that sylD1(f,d,x,y)[[1]] can be viewed as an alternate notation for a polynomial of degree d.

invSyl[S_,d_,x_,y_]:=FromCoefficientRules[Table[pExps[d][[i]]->S[[i]],{i,Length[S]}],{x,y}];

f=1+x^2-x y+y^3;
S=sylD1[f,3,x,y][[1]]

f1=invSyl[S,3,x,y]

We have the following identity: given polynomials f, g of degree, d1, d2 multiplication equation:

sylD1[f,d1,x,y].sylD1[g,d1+d2,x,y]=sylD1[f g,d1+d2,x,y]

g=2-x+x^2-3x y;
sylD1[f,3,x,y].sylD1[g,5,x,y]

Sfg=sylD1[f*g,5,x,y]

invSyl[Sfg[[1]],5,x,y]

h=Expand[f*g]

We can use the identity to perform numerical division. In general, one cannot divide one bivariate polynomial by another; the division algorithm does not directly apply to polynomials of several variables. However, if we know that some polynomial (say g) divides numerical polynomial h, then we can find another factor f using the multiplication equation above.

P=PseudoInverse[N[sylD1[g,5,x,y]]];
S=Chop[sylD1[h,5,x,y].P];
invSyl[S[[1]],3,x,y]

h2=-3.329102349826519` x-2.383334088473527` (x^2)-0.42613781708474535` (x^3)+2.2174484310530316` y-3.1536576003667216` x y-1.3649471149541397` (x^2) y-2.616665911526473` (y^2)-1.2169484594921673` x (y^2)-0.19250584655083047` (y^3)

ContourPlot[h2==0,{x,-8,2},{y,-2,8}]

Visually, it appears that this cubic is the union of a line and a quadric.

sol2={x,y}/.NSolve[{y^2-64,h2},{x,y},Reals]

g2= line[sol2[[1]],sol2[[4]],x,y]

P=PseudoInverse[N[sylD1[g2,3,x,y]]];
S=Chop[sylD1[h2,3,x,y].P];
f2=invSyl[S[[1]],3,x,y]

ContourPlot[{h2==0,f2==0},{x,-8,2},{y,-2,8}]

nDivide[h_,g_,x_,y_,tol_]:=Module[{d1,d2,P,S,f},
d1=tDeg[g,x,y];
d2=tDeg[h,x,y];
If [d1>d2,Print["Does Not Divide"];Return[Fail]];
P=PseudoInverse[N[sylD1[g,d2,x,y]],Tolerance->tol];S=Chop[sylD1[h,d2,x,y].P];f=invSyl[S[[1]],d2-d1,x,y];
If[Norm[Flatten[sylD1[f*g-h,d2,x,y]]]>d2*tol,Print[Style["Warning: g Does not divide h at this tolerance",Orange]];Return[Fail]];
f];

Our main goal of this subsection is our numerical greatest common divisor (GCD) algorithm. Suppose we have bivariate polynomials f and h with GCD g. Of course, g is only defined up to scalar multiple. If f and h have no common divisor of degree greater than 0, then g is a constant. Otherwise, we have bivariate polynomials u and v with f=g u and h=g v. Multiplying f by v and h by u givesf v=g u v=h u or f v-h u=0. Using the previous multiplication idea, some linear combination of the rows of certain large enough Sylvester matrices will be zero; hence there will be a nontrivial vector in the left null space. Since h=v and f=u satisfy the equation f v-h u=0 when the degree is tDeg(f)+tDeg(h), we need go no further than that. But if there is a nontrivial GCD, it will appear at a lower degree than that, so the following finitely terminating procedure will give the numerical GCD. Note that the simple While loop does the main work in this procedure.  We finish by noting that since u,v will be minimal then g=(f/u)=(h/v)is the GCD.  Because bivariate complex polynomials have the unique factorization property it actually follows that there will be an equation f v=h u of degree less than m+n if and only if f,h have a common factor.

Numerical GCDs ff and hh are bivariate polynomials in x, y. tol is the tolerance. The default is 10*^-12, but the larger the tolerance, the more likely it is to get a GCD.

nGCD[ff_,hh_,x_,y_,tol_]:=Module[{S,Sk,U,V,u,v,su,d1,d2,dv,f,h,n,m,k,W,lnull},
lnull[M_]:=Dimensions[M][[1]]-Length[SingularValueList[N[M],Tolerance->tol]];
d1=tDeg[ff,x,y];
d2=tDeg[hh,x,y];
If[d1≥d2,f=ff; h=hh;m=d1;n=d2,
f=hh;h=ff;m=d2;n=d1];
k=m-1;
Sk={{0}};
While[k<m+n && lnull[Sk]==0,
k++;                          (* Main loop *)
Sk=sylvesterMatrix[f,h,k,x,y]];
{U,S,V}=SingularValueDecomposition[N[Sk],Tolerance->tol];
W=Flatten[Take[U,All,-1]];
dv=k-m;
If[dv≤0, Return[h]];
su=Chop[Take[W,Binomial[dv+2,2]],tol];
u=invSyl[su,dv,x,y];
nDivide[h,u,x,y,tol]]

g=randomRealPolynomial[2,x,y]

h=Expand[g*randomRealPolynomial[3,x,y]]

f=Expand[g*randomRealPolynomial[4,x,y]]

gout=nGCD[f,h,x,y,1.*^-10]

Expand[g-gout/gout[[1]]*g[[1]]]

This algorithm was developed in 2003 independently by E. Kaltofen and me. My colleague Z. Zeng added a preprocessing stage to make the algorithm faster and a postprocessing stage to make it more accurate. His algorithm (Zeng and Dayton, 2004) is still used as a benchmark to test other GCD algorithms. For our purposes, the present version is fast and accurate enough.

A relevant use of the GCD algorithm in this book is to test the hypotheses of Bézout’s theorem, especially if the solver detects infinitely many solutions. Another application is to test to see if a bivariate polynomial is square-free and, if not, to give a square-free polynomial defining the same curve. Here, the theory says that if polynomial f is not square-free, then almost any partial derivative h=a D[f,x]+b D[f,y] will share a common divisor with f. Dividing f by the greatest common divisor will give a square-free polynomial defining the same curve.

f4=-39.72740279218701` (x^2)-6.576432679235679` (x^4)-0.27216394569483804` (x^6)+20.230248062255377` y+52.41558252827816` (x^2) y+4.199817230888195` (x^4) y-24.986060223213187` (y^2)-23.79488900809628` (x^2) (y^2)-0.5443278913896761` (x^4) (y^2)+11.063866289626182` (y^3)+4.338410336604424` (x^2) (y^3)-2.0680831293809843` (y^4)-0.27216394569483804` (x^2) (y^4)+0.13859310571622874` (y^5);

ContourPlot[f4==0,{x,-3,3},{y,-1,6}]

This fifth-degree curve looks like a parabola with some fuzz. Since ContourPlot works by changing the sign of f, any algebraic component that is a square will generally not show up but could cause some noise.

h4=Expand[{D[f4,x],D[f4,y]}.RandomReal[{-1,1},2]]

g4=nGCD[f4,h4,x,y,1.*^-10]

ff=nDivide[f4,g4,x,y,1.*^-12]

RegionPlot[ff>0,{x,-3,3},{y,-1,6}]

fp=nDivide[f4,g4^2,x,y,1.*^-12]

ContourPlot[{fp==0,g4==0},{x,-3,3},{y,-1,6}]

This solves our problem. However, in the algorithm, the last step was a division, but we found ff by another division, which just returned what we had originally, so it is better to have the square-free reduction as a separate procedure and save two divisions, as well as saving some clerical steps (because we know degree f is greater than degree h).

sqFree[f_,x_,y_,tol_]:=Module[{S,Sk,U,V,v,su,dv,h,n,m,k,W},
h=Expand[{D[f,x],D[f,y]}.RandomReal[{-1,1},2]];
lnull[M_]:=Dimensions[M][[1]]-Length[SingularValueList[N[M],Tolerance->tol]];
m=tDeg[f,x,y];
n=tDeg[h,x,y];
k=m-1;
Sk={{0}};
While[k<m+n && lnull[Sk]==0,
k++;
Sk=sylvesterMatrix[f,h,k,x,y]];
{U,S,V}=SingularValueDecomposition[N[Sk],Tolerance->tol];
W=Flatten[Take[U,All,-1]];
dv=k-n;
su=Chop[Take[W,-Binomial[dv+2,2]],tol];
invSyl[su,dv,x,y]];

sqFree[f4,x,y,1.*^-12]

SquareFreeQ[f4,{x,y}]

f5=4 y+5 (y^2)+7 x (y^2)+5 (x^2) (y^2)-2 (x^4) (y^2)+5 (y^3)+8 x (y^3);

ContourPlot[f5==0, {x,-4,4},{y,-4,4}]

Factor[f5]

The component y=0 only partly appears in the plot—it disappears, as in a squared component. Running sqFree with a tight tolerance returns the curve as square-free, but with a very loose tolerance, we get various answers because the algorithm uses randomization.

g=sqFree[f,x,y,.02]

ContourPlot[{g==0, f/y==0},{x,-4,4},{y,-4,4}]

So the almost-square component of f is removed, giving close to the second exact component of f.

A.4 Macaulay Matrices

Sylvester matrices are global objects because we use the whole defining polynomials, but Macaulay matrices are based at a point. Further, we will see that the degree of the polynomials is not used in defining Macaulay matrices (which allows us to extend the notion of intersection multiplicity to curves defined by analytic functions, as well as polynomial functions).

There are several strategies for calculating Macaulay matrices, which tend to be large. Fortunately, we will only work in two dimensions with polynomials of modest degree (where the number of entries in these matrices run in the thousands), but the concept can be expanded to any number of variables. The size of Macaulay matrices grows very quickly as the degree and number of variables increase, and entries are measured in millions. Thus in choosing methods, considering computer time and memory is important.

In choosing methods, I compared two strategies. The first, copying the idea of constructing Sylvester matrices, uses the fact that polynomials are given in terms of coefficients. We can use the fAssoc function in A.3 to get an association directly from CoefficientRules after using a change of basis calculation

fp=Expand[f/.Thread[{x,y}->{x,y}+p].

To calculate a single Macaulay matrix, this is twice as fast and adds about 60% less memory to the session.

However (as in the previous GCD algorithm), in our main intersection multiplicity algorithm that follows, we will again be increasing the degree until we reach a stopping point, but the Macaulay matrix of one degree larger has the previous Macaulay matrix as a submatrix, and one can possibly save time and memory by not recalculating these entries (although sometimes memory lookup can be more time consuming than recalculating). Thus the algorithm below (given the calculation of the previous degree Macaulay matrix) is 25% faster and adds about the same amount of new memory as the Sylvester algorithm. Although we won't need this, it will also work with analytics as well as polynomial functions.

We start with the differential operator difOp, which would calculate the coefficient of (x^i)(y^j) of the transformed function fp mentioned above:

difOp[f_,pt_,i_,j_]:=N[(1/Factorial[i]/Factorial[j]D[f,{x,i},{y,j}])/.Thread[{x,y}->pt]]

We assume that we will only calculate Macaulay matrices at common points of a pair of curves, so the functions (which can now be analytic as well as polynomials) will be approximately zero at the point pt. The recursive function dAss will be treated as a global variable, so all associations calculated will be stored for the remainder of the session.

dAss[f_,pt_,0]=<|{0,0}->0|>;
dAss[f_,pt_,n_]:=dAss[f,pt,n]=AssociateTo[dAss[f,pt,n-1],Table[{i,n-i}->difOp[f,pt,i,n-i],{i,0,n}]];

macaulayD1[f_,p_,n_]:=Module[{M},
Table[Values[shiftFA[dAss[f,p,n-Total[q]],q,n]],{q,pExps[n-1]}]]

macaulayMatrix[f_,g_,p_,n_]:=Join[macaulayD1[f,p,n],macaulayD1[g,p,n]]

f6=-0.007570349710845024` (x^2)+0.0038550230270411076` y;
g6=-72.4414468536883`-5.995940647212547` (x^2)+44.73564407546775` y-5.995940647212419` (y^2);
pt6={-1.3304175490338985`,3.475880965863065`}

M6=macaulayMatrix[f6,g6,pt6,3];
Dimensions[M6]

M6//MatrixForm

SingularValueList[M6]

{U,S,V}=SingularValueDecomposition[M6];
NS=Take[V,All,-2];
Chop[NS]//MatrixForm

diff1=NS[[1,2]]*difOp[#,pt6,0,0]&;
diff2=NS[[2,1]]*difOp[#,pt6,1,0]+NS[[3,1]]difOp[#,pt6,0,1]&;
{diff1[f6],diff1[g6]}

{diff2[f6],diff2[g6]}

We have two differential operators that numerically kill (tolerance (10^(-12))) both f and g. Leaving out the technical details, this is the numerical analyst’s definition of intersection multiplicity 2. (See the introduction to [Dayton, Li, and Zeng, 2011] for the details.)

We can now describe the intersection multiplicity algebraically according to F. S. Macaulay (1916). Given two curves f=0,g=0 and a point p on both curves, one computes the Macaulay matrices macaulayMatrix[f,g,p,i] of order i=1,2,3... calculating via SingularValueList the dimension n of the null space (that is, the number of columns minus the number of singular values larger than your tolerance). The default is (10^(-12)), but this can be changed. As i increases, this sequence, known as the Hilbert function, will be nondecreasing. When the same number appears twice, it will (theoretically) remain constant. This is the intersection multiplicity. The following procedure gives the multiplicity, optionally printing out the Hilbert function. It will assume the equations of the curves involve variables x and y, and no other variables appear. By default (to avoid an infinite loop), the procedure will stop at multiplicity maxm=20, but this is easy to change.

If the two equations have a nontrivial common divisor and the point p is on the common curve, then the multiplicity is infinity (that is, the procedure will not terminate). Intersection multiplicity of curves f,g are in variables x,y only at point p with tolerance tol.

Options[intersectionMultiplicity]={HF->False};
intersectionMultiplicity[f_,g_,p_,tol_,OptionsPattern[]]:=Module[{M,c,r,m,n,del, maxm,i,L,x,y},
maxm=20;
If[Abs[f/.Thread[{x,y}->p]]>tol,Print["p not zero of f"];Abort[]];
If[Abs[g/.Thread[{x,y}->p]]>tol,Print["p not zero of g"];Abort[]];
del=1;
m=1;
L={1};
For[i=1,i<maxm && del>0,i++,
M=macaulayMatrix[f,g,p,i];
c=Dimensions[M][[2]];
r=Length[Select[SingularValueList[M],#>tol&]];
n=c-r;
del=n-m;
m=n;
AppendTo[L,m]];
If[m==maxm,Print["Multiplicity Greater or equal to ",maxm]];If[OptionValue[HF],Print["depth = ",i-2," HF = ",L]];
m]

Some examples follow.

dTol=1.*^-9;
intersectionMultiplicity[x,y,{0,0},dTol]

intersectionMultiplicity[y-x^2, x-y,{1,1},dTol]

intersectionMultiplicity[y,y-x^2,{0,0},dTol]

intersectionMultiplicity[x-y+x^7,x-y + y^7,{0,0},dTol, HF->On]

intersectionMultiplicity[x^3, x^2 y+y^4,{0,0},dTol,HF->On]

This algorithm was originally proposed by F. S. Macaulay in 1916 but was largely forgotten because it required finding the rank of large matrices. My colleague Z. Zeng and I implemented his algorithm numerically in 2004 (a discussion on the topic can be read in [Dayton, Li, and Zeng, 2011]). A problem that could arise is when an intersection has high multiplicity, functions such as NSolve may not be able to calculate the point accurately, and then intersectionMultiplicity may not give the correct multiplicity. Fortunately, at approximately the same time, our neighbors A. Leykin, J. Verschelde, and A. Zhao (2006) proposed a deflation algorithm that would replace the system in the non-quadratically converging Newton’s method with a quadratically converging system. These two algorithms are closely related and, ideally, should be used together. Unfortunately, in the case that the multiple intersection point is not known exactly, the deflation algorithm may be too complicated for normal use.

f8=-4+12 x-27 (x^2)+27 (x^3)-3 y;
g8=5+3 x+33 y+54 (y^2)+27 (y^3);
p8={x,y}/.NSolve[{f8,g8},{x,y},Reals][[1]]

intersectionMultiplicity[f8,g8,p8,1.*^-9]

J8=Det[D[{f8,g8},{{x,y}}]]

p8a=gaussNewton3[f8,g8,J8,p8,x,y]

intersectionMultiplicity[f8,g8,p8a,1.*^-9]

p8e=N[{1/3,-2/3}]
intersectionMultiplicity[f8,g8,p8e, 1.*^-9, HF->On]

intersectionMultiplicity[f8,g8,p8,.02]

A.5 Finding Algebraic Equations of Curves Parameterized by Polynomials

In Chapter 1, we mentioned that curves parameterized by polynomial or rational parameterizations are contained in algebraic curves. In Chapter 2, we suggested finding these parameterizations by choosing an appropriate number of random points on the curve and using interpolation to find the algebraic curve.

Here, we show that there is another approach for polynomial parameterizations that does not require finding points. This is especially advantageous if the polynomials are exact (integer coefficients), as the algebraic curve is exact. Moreover, we can explain why it is always possible to find an algebraic equation of a curve containing the parameterized curve. Finally, we will sketch (but not give details) on how to extend this to rational parameterizations.

The main tool is a matrix I will call the p2aMatrix, which is somewhat analogous to the Sylvester matrix. For convenience, we will index the rows by monomials of total degree d or less and the columns by powers of t. Assuming the parameterized curve is given by F={u,v} where u, v are polynomials in t, then the entry in the row (x^i)(y^j) and column (t^k) is the coefficient of (t^k) in the expansion of (a^i)(b^j). Since the row indices have total degree d or less, (a^i)(b^j) will have as a polynomial in t degree no larger than m d, where m is the larger of the degrees of a, b. Thus the p2aMatrix will have dimensions (((d+2)(d+1))/2)×(d m+1). To get a result, we will need d≥m.

F1={t^2,2t-1};

The maximum degree is m=2, so we try d=2. Then upon substitution of the parameters into the monomials, and if the constant remains 1, x becomes (t^2); y becomes 2 t-1; (x^2) becomes (t^4); xy becomes (t^2)(2t-1)=2(t^3)-(t^2); and (y^2) becomes (((2t-1))^2)=4(t^2)-4t+1.

Thus the p2aMatrix is 6×5, and by elementary linear algebra, the rows must be dependent, that is, some linear combination of the rows gives the zero row. By inspection, we see that if we add row 1, -4 times row x, twice row y, and row (y^2), we get the zero row. In other words, we could write 1-4x+2y+(y^2)=0.

Expand[(1-4x+2y+y^2)/.Thread[{x,y}->F1]]

That means every point on the parameterized curve lies in the curve 1-4x+2y+(y^2)=0.

By this approach, we can find a curve that contains our parameterized curve. In the next example, the curve is a parabola and every point in the parabola can be written in the form F(t) for some t. In general, the parameterized curve may lie in some topological component of the algebraic curve and need not give the entire curve.

Show[ParametricPlot[F1,{t,-3,3},Axes->False,Frame->True, PlotStyle->Cyan],    ContourPlot[1-4x+2y+y^2==0,{x,-1,9},{y,-7,5}, ContourStyle->Dashed],ImageSize->Medium]

p2aMatrix[F_,d_,t_,x_,y_]:=Module[{m,n,M,c,ca,ra},
m=Max[Exponent[F[[1]],t],Exponent[F[[2]],t]];
n=m*d;
M={};
Do[c=CoefficientRules[Expand[(x^p[[1]] y^p[[2]])/.Thread[{x,y}->F]],t];
ca=Association[c];
ra=Association[Table[If[MissingQ[ca[{k}]],{k}->0,{k}->ca[{k}]],{k,0,n}]];
AppendTo[M,Values[ra]],{p,pExps[d]}];
M];

p2aMatrix[F1, 2,t,x,y]//MatrixForm

The theory is that since m is fixed by the parameterization, the number of rows is a quadratic function of d, but the number of columns is only a linear function of d. Thus for some d, we must have more rows than columns. There will be a linear combination of the rows adding to zero, and our parameterized curve will be contained in the corresponding algebraic curve. We will see that we might get such a linear combination even before this happens. Using the smallest value of d will make it likely that the algebraic curve will not be reducible with algebraic components (which are not needed).

It is common sense that our algebraic equation should have degree s=m since the parameterization has a (t^( m)) in it. We will see that this isn’t actually true, but counterexamples are unlikely (and it will be apparent if we have one). We calculate p2aMatrix for d=m and look for a nonzero vector u in the left null space. How we do it will depend on whether we want an exact or numerical answer. In either case, we can use the following procedure to find the equation.

p2aEq[u_,x_,y_]:=Module[{P,d},
If[Length[Dimensions[u]]>1,Print["u must be a vector"];Abort[]];
d=Switch[Length[u],3,1,6,2,10,3,15,4,21,5,28,6,36,7,45,8,55,9];
If[Length[d]>1,Abort[]];
P=pExps[d];
FromCoefficientRules[Table[P[[k]]->u[[k]],{k,(d+2)(d+1)/2}],{x,y}]];

If our parametric function has integer coefficients for both component polynomials, then p2aMatrix will be an integer matrix and Mathematica will treat it as such. Calculating the left null space of PM=p2aMatrix(F,d,t,x,y) by NullSpace[Transpose[PM]] will result in (at worst) rational numbers, which can be converted to integers by multiplying by the greatest common denominator.

c5={ChebyshevT[4,t],ChebyshevT[5,t]}

PM=p2aMatrix[c5,5,t,x,y];

NS=NullSpace[Transpose[PM]]

f=p2aEq[NS[[1]],x,y]

Expand[f/.Thread[{x,y}->c5]]

Show[ParametricPlot[c5,{t,-2,2},Axes->False, Frame->True,PlotRange->2,PlotStyle->Cyan],ContourPlot[f==0,{x,-2,2},{y,-2,2}, ContourStyle->Dashed]]

G1={t^2-1,t^2+1};
PG1=p2aMatrix[G1,2,t,x,y];

NS=NullSpace[Transpose[PG1]]

g1a=p2aEq[NS[[1]],x,y]
g1b=p2aEq[NS[[2]],x,y]
g1c=p2aEq[NS[[3]],x,y]

{ContourPlot[g1a==0,{x,-4,4},{y,-4,4}],ContourPlot[g1b==0,{x,-4,4},{y,-4,4}],ContourPlot[g1c==0,{x,-4,4},{y,-4,4}], ParametricPlot[G1,{t,-2,2},Axes->False,Frame->True, PlotRange->4]}

In this case, the parametric curve is half the line given by g1c. The other two curves have an extra algebraic component, as well as the extra half line. The curve has degree 1 while m=2.

PG1a=p2aMatrix[G1,1,t,x,y];
NS1=NullSpace[Transpose[PG1a]]

p2aEq[NS1[[1]],x,y]

F=N[{CubeRoot[2]t^2,t^3}]

PF=p2aMatrix[F,3,t,x,y];
Dimensions[PF]

{u,s,v}=SingularValueDecomposition[PF];
Diagonal[s]

ns=Flatten[Take[u,All,-1]]

f=p2aEq[ns,x,y]

Expand[f/Coefficient[f,x^3]]

The previous examples show that while integer polynomial parameterized curves are contained in exact curves, exact curves most often do not have integer polynomial parameterizations but may have numerical parameterizations.

Finally, this method can be extended to deal with rational parameterizations. We can assume the two component functions have the same denominator, i.e. F={(a/c), (b/c)} for polynomials a,b,c in t. Then by homogenizing, we get a polynomial parameterization H={a,b,c}. Extending the procedures p2aMatrix, p2aEq to three variables and adding a random second equation (Subscript[r, 1])x+(Subscript[r, 2])y+(Subscript[r, 3]) z=1, we can find a homogeneous equation h[x,y,z]=0 whose projective solution contains the image of H. Setting z=1, for example, gives us an affine equation.

A.6 Factoring Highly Singular Curves

In Chapter 5, we gave the famous classical result that an irreducible curve of degree d has, at most, (((d-1)(d-2))/2) singularities. Rather than the usual indirect proof, we sketched the direct proof of the contrapositive: a curve of degree d with more than (((d-1)(d-2))/2) singularities is reducible. I promised that we could actually make this into an algorithm implemented in Mathematica. We do this here.

Assume that at least (((d-1)(d-2))/2)+1 is affine (possibly complex).

singularFactor[f_,x_,y_,tol_]:=Module[{h,csp,d,sp,M,M1,cr,c0,rhs,rr,XX,h1},
csp=Chop[complexSingularPoints[f,x,y]];
d=tDeg[f,x,y];
If[d<3,Print["degree must be at least 3"];Abort[]];
If[Length[csp]<((d-1)(d-2)/2+1),Print["Not enough affine singular points"];Abort[]];
csp=RandomSample[csp,(d-1)(d-2)/2+1];
sp=If[d==3,csp,Join[csp,realRegularPoints[f,x,y,d-3]]];
XX=pExps[d-2];
M=Table[Which[e=={0,0},1,e[[1]]==0,p[[2]]^e[[2]],e[[2]]==0,p[[1]]^e[[1]],True,p[[1]]^e[[1]]*p[[2]]^e[[2]]],{p,sp},{e,XX}];
rr=Table[0,{Length[XX]}];
rr[[-1]]=1;
M1=Append[M,RandomReal[{-1,1},Length[XX]]];
rhs=Table[0,{Length[XX]}];
rhs[[-1]]=1;
cr=LeastSquares[M1,rhs];
c0=Min[Select[Abs[cr],#>.1&]];
h=FromCoefficientRules[Table[XX[[i]]->cr[[i]],{i,Length[XX]}],{x,y}];
h=ComplexExpand[Re[h]];
h1=nGCD[f,h,x,y,tol];
If[tDeg[h1,x,y]==0,Print[Style["Warning: GCD not found, Returning h",Orange]];Return[{h,1}]];
Chop[{h1,nDivide[f,h1,x,y,tol]},10*tol]];

g7=18-24 x+12 (x^2)-8 (x^3)+2 (x^4)-24 y+32 x y-8 (x^2) y+9 (y^2)-8 x (y^2)+3 (x^2) (y^2)-4 (y^3)+(y^4);

complexSingularPoints[g7,x,y]

{h,k}=singularFactor[g7,x,y,1.*^-8]

ContourPlot[{h==0,k==0}, {x,-2,4},{y,-2,4}]

f7=(x^2)+5 (x^3)+8 (x^4)+4 (x^5)-(y^2)-4 x (y^2)-4 (x^2) (y^2)+(x^2) (y^3)+(x^3) (y^3)-(y^5);
rsp=Chop[realSingularPoints[f7,x,y]]

ContourPlot[f7==0,{x,-2,2},{y,-3,3}, Epilog->{PointSize[Medium],Red,Point[rsp]}]

Length[complexSingularPoints[f7,x,y]]

Since (((6-1)(6-2))/2)+1=11, f should factor. However, the curve through fourteen points is of degree 4, and there is no factor of degree 4 or complementary degree 2.

{h,k}=singularFactor[f7,x,y,1.*^-17]

ContourPlot[{f7==0,h==0},{x,-3,3},{y,-5,5}]

h contains an extra line. This is why singularFactor needs the nGCD function.

{h,k}=Chop[singularFactor[f7,x,y,1.*^-7],1.*^-7]

ContourPlot[{f7==0,h==0},{x,-3,3},{y,-5,5}]

singularFactorInf[f_,x_,y_,tol_]:=Module[{d,fh,f1,z,A,h,k},
d=tDeg[f,x,y];
A=iTransform[RandomReal[{-1,1},3].{x,y,1},x,y];
f1=FLT[f,A,x,y];
{h,k}=singularFactor[f1,x,y,tol];
Chop[{FLT[h,Inverse[A],x,y],FLT[k,Inverse[A],x,y]},10*tol]
];

g9 =5-2 x+2 y-6 x y+2 (x^2) y-2 x (y^2)+(x^2) (y^2);

Length[complexSingularPoints[g9,x,y]]

{h,k}=singularFactorInf[g9,x,y,1.*^-7]

ContourPlot[{h==0,k==0},{x,-5,5},{y,-5,5}]

This is not a general factoring algorithm as there are many reducible curves of degree d>3 with fewer than (((d-1)(d-2))/2) singular points.

In fact, it seems that the only examples of reducible curves with more than this many singular points are curves such that each factor has almost the maximum number of singular points. For example, reducible curves with factors of degree 1 and 2, singular curves of degree 3, and Chebyshev curves.

A.7   A note on Bézout’s Theorem

(Added January 2019)

One omission in my book is a proof of Bézout’s theorem.  At the time  all existing proofs that I knew of used concepts that I have avoided in the book so this omission was intentional.  A problem is that my version of intersection multiplicity, which is included in the book specifically for the purposes of the statement of this theorem, is incompatible with these proofs, in particular to the ones by Abyhankar  [1] and Fischer [8] that I referenced in the book.

Recently I became aware of a paper by Telen, Mourrain and Van Barel in [Siam J. Matrix Anal Appl, vol 39 no3, pp 1421-1447,(2018)] giving a solution method for two variable polynomial systems which is related to material in Appendix 1, in particular my definition of  intersection multiplicity can apply.  Although what follows is not a  proof of the theorem, it does illustrate the connection between the solutions, including complex, multiple and infinite, and the nullspace of a Sylvester matrix.

My argument involves computing the nullspace of the Sylvester matrix of order m+n of relatively prime polynomials of degree m,n in two different ways and comparing.  If you are using the notebook version please initialize global functions and this notebook.  There will be some new functions defined here but they will not be added to the global functions page since they are not elsewhere in the book.

Clear[f1,f2,f3,f4,g1,g2,g3,g4]

The nullity of the Sylvester Matrix

Let f,g be the two bivariate polynomials of degrees m,n respectively.  Our discussion in section A.3 concludes that if GCD[f,g]=1, that is, if  f,g have no common component (factor), then the Sylvester Matrix of order  m+n-1 or less will have independent rows, that is the row rank will be the number of rows.  But since f g and g  f will correspond to two different rows of the Sylvester matrix of order m+n the rank of the rowspace will be only the number of rows minus 1.

The nullity, that is the vector space dimension of the null space, will be the number of columns minus the row rank.  The number of columns of the Sylvester matrix of order m+n will be the number of monomials of degree less than or equal to m+n, calculated in Chapter 2.

cols=Expand[(n+m+2)(m+n+1)/2]

The row rank from the above will be the number of monomials of degree ≤ m plus the number of monomials of degree ≤n  minus 1.  Thus it will be

rrank=Expand[((n+2)(n+1)+(m+2)(m+1))/2-1]

Subtracting

cols-rrank

This, and the unique factorization property actually prove the following Theorem :

Little Bézout Theorem :  Let f,g be bivariate polynomials, of degrees m,n respectively, with no common factor.  Then the null space of the Sylvester matrix of order  m+n of  f,g has dimension m n.

Below we will identify this null space with the intersection points of the two curves, i.e. the solution set of the polynomial system {f,g} in the complex projective plane.  We will work backwards from the solutions so this will not be a proof of Bézout’s theorem but an explanation of what is happening.

Case of  simple affine solutions

Each affine point will correspond to an evaluation vector, or, e-vector.  We define two simple functions, these will not be listed in Appendix 2 since they are only used here.

pMons[d_,x_,y_]:=Table[x^p[[1]] y^p[[2]],{p,pExps[d]}]

Example:

pMons[3,x,y]

Our evaluation vectors are now defined for a point p in (ℂ^2)and a degree d by evaluating each monomial at the same given point.

eVec[p_,d_]:=pMons[d,x,y]/.Thread[{x,y}->p];

Example :

e1=eVec[{3,2},4]

For simplicity we will use an example from the article TMVB cited above since it has only integer solutions.

f1=7+3x-6y-4x^2+2x y+5y^2;
g1=-1-3x+14y-2x^2+2x y-3y^2;

sol1={x,y}/.NSolve[{f1,g1}]

S1=sylvesterMatrix[f1,g1,4,x,y];
NS1=Transpose[NullSpace[S1]];
NS1//MatrixForm

Thus the Little Bézout theorem correctly tells us there are 4 columns.  Calculating evaluation sequences for the other solutions

e2=eVec[{2,1},4]
e3=eVec[{-1,0},4]
e4=eVec[{-2,3},4]

S1.e1
S1.e2
S1.e3
S1.e4

We get 0 in each case!  Note one technicality, since Mathematica does not distinguish lists as being rows or columns it writes both rows and columns when given separately as rows.  When multiplying it switches rows to columns by the context.  So here the Sylvester matrix row is treated as a row while the evaluation sequence row is treated as a column.

So these vectors as columns are in the null space.  We already see that e3 is the last column of the nullspace as calculated by Mathematica’s NullSpace function.  But if we look at

NS1a=Transpose[{e1,e2,e3,e4}];
NS1a//MatrixForm
MatrixRank[NS1a]
S1.NS1a

we see that NS1a is a matrix of rank 4 with columns in the null space of S1 and thus the columns give an alternate basis for the null space of S1.

It is not hard to prove in general that degree d evaluation vectors of solution points of a system {f,g} lie in the null space of the Sylvester matrix of f,g of order d.   It is harder to prove directly that these e-vectors of distinct solution points are independent but will not be proven here.

The case of affine multiple points

If there are affine multiple points, we now show to find additional vectors in the nullspace of the Sylvester matrix.  Consider the following Example 2

f2=-4 x+4 (x^2)+(y^2);
g2=x^2+y^2-1;

The point {1,0} is easily seen to be on both curves.

intersectionMultiplicity[f2,g2,{1,0},1.*^-12]

From our Dayton-Li-Zeng paper  [6] on multiplicity this tells us there will be two local dual differentials at that point, or in other words two independent vectors in the nullspace of the Macaulay matrix at that point.  Again we will  take the order of the Macaulay matrix to be the sum of the degrees

M4a=macaulayMatrix[f2,g2,{1,0},4];
 NM4a=Chop[NullSpace[M4a],1.*^-14]

However we are interested in the Sylvester matrix which is not dependent on a point.  Fortunately there is a transformation which takes null vectors from a Macaulay matrix to the Sylvester matrix.  I won’t try to fully explain this here, but think of expanding a Taylor series about a point other than zero to get a series based at zero and the binomial coefficients play a large role.  This generalizes to the multivariable case and will be discussed in my Space Curve book.  In this function q is the point and n is the order of the Macaulay, Sylvester matrices.  This is strictly a matrix operation so no variables are involved

c2zMat[q_,n_]:=Module[{m,Tn,ss,bi,bj,r1,C,s,pow},
pow[a_,m_]:=If[m<=0,1,a^m];
s=Length[q];
Tn=pExps[n];
ss=Length[Tn];
ss=Length[Tn];
C={};
Do[bj=Tn[[j]];
C=Append[C,Table[Product[Binomial[Tn[[i]][[k]],bj[[k]]]*pow[q[[k]],(Tn[[i]][[k]]-bj[[k]])],{k,s}],
{i,ss}]],
{j,ss}];
Transpose[C]]

Dimensions[c2zMat[{1,0},4]]

NM4a0=c2zMat[{1,0},4].Transpose[NM4a];
NM4a0//MatrixForm

S4a=sylvesterMatrix[f2,g2,4,x,y];
S4a.NM4a0

So NM5a0 gives 2 independent vectors in the nullspace of S5.  Note also that

u1=eVec[{1,0},4]

is one of these vectors.

Our calculation of the dimension of the nullspace in the previous section had nothing to do with the multiplicity of the zeros.  There are 4 independent vectors since f2 and g2 are of degree 2   We check anyway.

MatrixRank[NullSpace[S4a]]

Solving

sol2={x,y}/.NSolve[{f2,g2}]

Then adding the evaluation vectors for these first two points gives

K=Join[NM4a0,Transpose[{eVec[sol2[[1]],4],eVec[sol2[[2]],4]}],2];
Dimensions[K]
MatrixRank[K]

Chop[S4a.K]

So we can conclude that K is equivalent to the null space of S5a.

As in the earlier section the independence generalizes.  Again we will not prove this but assure the reader that it is true.

The case of infinite points

If the rank of the matrix constructed from the affine points, simple and multiple, is not m n this suggests infinite points.  In this subsection we compute additional vectors in the nullspace from the infinite points.   The key is our calculation of infinite points in Chapter 3 of the book.  The infinite points of a curve come from solutions of the of maximal form of  the equation of the curve.  Then the infinite points of the intersection of two curves come from common solutions of the maximal forms of both equations.

I should mention that unless there is some geometrical reason why two curves should share infinite points, that infinite points are rare.  In particular if f,g are random curves it is very unlikely that there will be an infinite point in the intersection.  Thus our examples  are somewhat contrived.

These common forms appear, possibly after multiplication by other polynomials, in the upper degrees of the order m+n Sylvester matrix.    For simple, multiplicity 1, crossings one need only look at the highest order.  It is easily seen that the number of monomials in two variables of degree exactly m+n is m+n+1 so we find the nullspace of the matrix given by the last  m+n+1 columns of the Sylvester matrix of order m+n.  Unlike the affine case we handle the multiple case of infinite points together with the simple case but we may need to look in a lower order as well.

A simple example is the intersection of two rational function curves y=(1/(x-2)), y=((-1)/(x+2))

f3=Expand[y(x-2)-1]
g3=Expand[y(x+2)+1]

ContourPlot[{f3==0,g3==0},{x,-5,5},{y,-10,10},ImageSize->Small]

NSolve[{f3,g3}]

These are both quadratics so we expect 4 intersection points, hence 3 must be infinite.

infiniteRealPoints[f3,x,y]
infiniteRealPoints[g3,x,y]

Two are found, essentially {1,0,0}, {0,1,0}, but one must be multiple.  One can check using the method of Chapter 8 but from the plot it appears that {1,0,0} is the culprit, which is actually the case.  Calculating the Sylvester matrix and its nullspace gives

S3=sylvesterMatrix[f3,g3,4 ,x,y];
Dimensions[S3]
ns3=Transpose[NullSpace[S3]];
ns3//MatrixForm

with 4 columns as expected.  But finding the nullspace of the last 5 columns of S3, the degree 4 part, gives only

nsinf=Transpose[NullSpace[Take[S3,All,-5]]];
nsinf//MatrixForm

So we try including the degree 3 part which has 4 columns.

nsinf=Transpose[NullSpace[Take[S3,All,-9]]];
nsinf//MatrixForm

This has the desired rank 3 corresponding to 3 infinite points by multiplicity.  Since the actual Sylvester matrix S11 has 15 columns we prepend a 6×3 matrix of zeros

infpart= Join[Table[0, {i,6},{j,3}],nsinf];
Dimensions[infpart]

Finally we can add the evaluation vector corresponding to the affine solution {0, -.5} as a column

ev1=Partition[eVec[{0,-.5},4],1];
nullMat=Join[ev1,infpart,2];
nullMat//MatrixForm

We see

MatrixRank[nullMat]
S3.nullMat

so nullMat gives a basis for the nullspace of S11 once again explaining this nullspace in terms of the intersection points of the two curves.

Again I assert, but will not prove, that this example is typical.

One Final example

One situation where we can expect infinite intersection points is with gaussCurves of  monic polynomials of the same degree.

f4=gaussCurve[z^3-3z^2,x,y,z]
g4=gaussCurve[z^3+2z^2-z,x,y,z]

ContourPlot[{f4==0,g4==0},{x,-10,10},{y,-10,10}]

Normalize/@infiniteRealPoints[f4,x,y]
Normalize/@infiniteRealPoints[g4,x,y]

intersectionMultiplicity[f4,g4,{0,0},1.*^-12]

So here we have both common infinite points but also a multiple intersection point at the origin.  The affine intersection points are all complex, except for the origin.

sol={x,y}/.NSolve[{f4,g4}]

S6=sylvesterMatrix[f4,g4,6,x,y];

NS6=Transpose[NullSpace[S6]];
NS6//MatrixForm

Dimensions[S6]

For the multiple point at the origin we find the Macaulay matrix

M6=macaulayMatrix[f4,g4,{0,0},6];
MVecs=Chop[Transpose[NullSpace[M6]]]

There is no transformation c2zMat necessary because the multiple zero is {0,0}  and c2zMat[{0,0},6] is the identity.  For the 3 infinite points which must be simple

S6E=Take[S6,All,-7]

N6E=Transpose[NullSpace[S6E]];
N6E//MatrixForm

Ivecs=Join[Table[0,{i,21},{j,3}],N6E]

S6.Ivecs

EVects=Transpose[Table[eVec[sol[[i]],6],{i,4}]];

To see the full nullspace from the zeros you may need to expand this notebook to wider viewing width .

NS6A=Join[EVects,MVecs,Ivecs,2];
NS6A//MatrixForm

Chop[S6.NS6A]

MatrixRank[NS6A]

This may not look like the nullspace directly from the Sylvester matrix but remember that this nullspace is a complex vector space.  It is not too difficult to find a complex tranisition matrix between this two matrices which give alternate bases of the same complex space.

Conclusion

From the discussion above in all the examples we were able to completely describe the null space of the Sylvester matrix of order m+n in terms of the calculated solutions.   While this is not an actual proof of Bézout’s theorem this should give some insight as to why it works.  A complete proof would need to theoretically prove independence of these columns and also show existence of these intersections from the nullspace.  The cited paper of Telen, Mourrain and Van Barel above does show this existence for at least the affine solutions.  On the other hand, working backwards from Bézout’s theorem we see that our claim of the existence and independence of this basis for the nullspace of the Sylvester matrix is supported.

As a bonus, we see  that to calculate the intersection multiplicity we only need to calculate one Macaulay matrix, that of order m+n and find the dimension of its nullspace.  We propose a new algorithm

intersectionMultiplicity2[f_,g_,p_,x_,y_]:=Module[{m,n,M,c,svdl},
m=tDeg[f,x,y];
n=tDeg[g,x,y];
M=macaulayMatrix[f,g,p,n+m];
c=Dimensions[M][[2]];
svdl=SingularValueList[M,Tolerance->1.*^-12];
c-Length[svdl]]

{t1,i1}=Timing[intersectionMultiplicity2[f2,g2,{1,0},x,y]]
{t2,i2}=Timing[intersectionMultiplicity[f2,g2,{1,0},1.*^-10]]

Unsurprisingly, the new algorithm may not be faster when the multiplicity is small, but look what happens when the multiplicity (and depth) is large.

{t3,i3}=Timing[intersectionMultiplicity2[x-y+x^8,x-y+y^8,{0,0},x,y]]
{t4,i4}=Timing[intersectionMultiplicity[x-y+x^8,x-y+y^8,{0,0},1*^-10]]

I suspect that this bound m+n of the depth (see our paper [6]) of a multiple point of the solution of a 2×2 polynomial system generalizes to  N×N systems.  That is, if the degrees of the N equations are (Subscript[d, 1,])…, (Subscript[d, N]) then the bound on the depth is D=(Subscript[d, 1])+⋯+(Subscript[d, N]).  Unfortunately even this bound is too large to be of much help in calculating intersection multiplicity for N>2 so such a result would only be of theoretical interest.

Tech Notes